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Using the Martin-Siggia-Rose method, we study propagation of acoustic waves in strongly heteroge- 
neous media which are characterized by a broad distribution of the elastic constants. Gaussian- white 
distributed elastic constants, as well as those with long-range correlations with non-decaying power- 
law correlation functions, are considered. The study is motivated in part by a recent discovery 
that the elastic moduli of rock at large length scales may be characterized by long-range power-law 
correlation functions. Depending on the disorder, the renormalization group (RG) flows exhibit a 
transition to localized regime in any dimension. We have numerically checked the RG results using 
the transfer-matrix method and direct numerical simulations for one- and two-dimensional systems, 
respectively. 



Understanding how waves propagate in heterogeneous 
media is fundamental to such important problems as 
earthquakes, underground nuclear explosions, the mor- 
phology of oil and gas reservoirs, oceanography, and med- 
ical and materials sciences [1] . For example, seismic wave 
propagation and reflection are used to not only estimate 
the hydrocarbon content of a potential oil or gas field, but 
also to image structures located over a wide area, ranging 
from the Earth's near surface to the deeper crust and up- 
per mantle. The same essential concepts and techniques 
are used in such diverse fields as materials science and 
medicine. 

In condensed matter physics, a related problem, 
namely, the nature of electronic states in disordered ma- 
terials, has been studied for several decades and shown 
to depend strongly on the spatial dimensionality d of 
the materials [2]. It was rigorously shown that, for one- 
dimensional (ID) systems, even infinitesimally small dis- 
order is sufficient for localizing the wave function, irre- 
spective of the energy [3], and that the envelop of the 
wave function ip{r) decays exponentially at large dis- 
tances r from the domain's center, ip(r) ~ exp(— r/£), 
with £ being the localization length. The most important 
results for d > 1 follow from the scaling theory of local- 
ization [4] which predicts that, for d < 2, all electronic 
states are localized for any degree of disorder, while a 
transition to extended states - the metal-insulator tran- 
sition - occurs for d > 2 if disorder is sufficiently strong. 
The transition between the two states is characterized by 
divergence of the localization length, £ oc \W — W c \~ v , 
where W c is the critical value of the disorder intensity. 



Wegncr [5] derived a field-theoretic formulation for the lo- 
calization problem which, together with the scaling the- 
ory [4], predict a lower critical dimension, d c = 2, for 
the localization problem. These predictions have been 
confirmed by numerical simulations [6] . 

Wave characteristics of electrons suggest that the lo- 
calization phenomenon may occur in other wave propa- 
gation processes. For example, consider propagation of 
seismic waves in heterogeneous rock. In this case, the 
interference of the waves that have undergone multiple 
scattering, caused by the heterogeneities of the medium, 
may cause their localization. Unlike electrons, however, 
the classical waves do not interact with one another 
and, therefore, propagation of such waves in heteroge- 
neous media (such as porous rock) provides [7] an ideal 
model for studying the classical Anderson localization 
[8-12] in strongly disordered media. This is the focus 
of this Letter. We study localization of acoustic waves 
in strongly heterogeneous media, and formulate a field- 
theoretic method to investigate the problem in the media 
that are characterized by a broad distribution of the elas- 
tic constants. Localization of acoustic waves was previ- 
ously studied by several groups [13], although not in the 
strongly disordered media that we consider in this Letter. 
The system that we study is the continuum limit of an 
acoustic system with off-diagonal disorder. Our approach 
is based on the method first introduced by Martin, Siggia 
and Rose [14] for analyzing dynamical critical phenom- 
ena. We calculate the one-loop beta functions [8,14] for 
both spatially delta-correlated and power-law correlated 
disorder in the elastic constants, and show that in any 



case there is a disorder-induced transition from delocal- 
ized to localized states for any d. In addition to be inter- 
esting on its own, our study is motivated in part by the 
recent discovery [15] that the distribution of the elastic 
moduli of heterogeneous rock contains long-range corre- 
lations characterized by a nondecaying power-law corre- 
lation function. Baluni and Willemsen [16] studied prop- 
agation of acoustic waves in a ID layered system, which 
can be thought of as a simple model of rock (although 
their goal was not to study acoustic wave propagation 
in rock), and showed that the waves are localized. How- 
ever, they did not consider higher-dimensional systems, 
nor did they study the type of disorder that we consider 
in the present Letter. The possibility of wave localiza- 
tion in disordered media with long-range, nondecaying 
correlations has important practical implications which 
we will discuss briefly. However, our results are com- 
pletely general and apply to any material in which the 
local elastic constants are distributed broadly. We con- 
firm the analytical predictions for ID and 2D systems 
using numerical simulations. 

Wave propagation in a medium with a distribution of 
clastic constants is described by the following equation 
(for simplicity we consider the scalar wave equation): 



d 2 



V-[A(x)VV(x,i)]=0, 



(1) 



where ip(x, t) is the wave amplitude, and A(x) = 
C(x)/m is the ratio of the elastic stiffness C(x) and the 
mean density m of the medium. We then write A as, 



A(x) = A + ?y(x) , 



(2) 



where Ao = (A(x)). In this Letter we assume 77 (x) to 
be a Gaussian random process with a zero mean and the 
covariancc, 



(r ? (x)r ? (x'))= 2A'(|x - x'|) = 2D S d ( X - x') 
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in which Do and D p represent the strength of the dis- 
order due to the delta-correlated and power-law corre- 
lated parts of the disorder. Previously, Souillard and 
co-workers [17] studied wave propagation in disordered 
fractal media, which is characterized by a decaying power- 
law correlation function. Their study is not, however, di- 
rectly related to our work. Consider a wave component 
■0(x, tu) with angular frequency lo, which is obtained by 
taking the temporal Fourier transform of Eq. (1) which 
yields the following equation for propagation of a wave 
component in a disordered medium, 



Since ?/(x) is a Gaussian variable, we obtain a Martin- 
Siggia-Rosc effective action S e for the probability density 
functional of the wave function ip(x, uj), given by 

dxdx'[(i7Mx')(V 2 + -r-)^j(x) 

+i^(x')(V 2 + ^)^(x). 

+X*(x')(V 2 + ^) X (x)) ( 5(x-x') 
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Here, tpiipc), %j) b(x), \ and \* are t ne auxiliary and 
Grassmanian fields of the field-theoretic formulation, re- 
spectively. Two coupling constants, go = Dq/Xq, and, 
g p = Dp/\ 2 , appear in S e . Thus, we carry out a 
renormalization group (RG) analysis in the critical limit, 
w 2 /Ao — > 0, to derive, to one-loop order, the beta func- 
tions [7,12] that govern the two couplings under the RG 
transformation. The results arc given by, 



dgo 
dlnl 

dg P 
dlnl 



= -dg 



10g 2 + 20.9o.9p 



(6) 



= (2p - d)g p + I2g g p + 16g p , (7) 



where I > 1 is the re-scaling parameter, and go and g p 
are given by, 



9 P = kd 




(8) 



(9) 



with kd = Sd/ (27r d ), and Sd being the surface area of the 
d— dimensional unit sphere. Examining the RG flows, 
Eqs. (6) and (7), reveals that, depending on p, there are 
two distinct regimes: 

(i) For < p < d/2 there are three sets of fixed points: 
The trivial Gaussian fixed point (c/q = g* = 0) which is 
stable, and two non-trivial fixed points and eigendirec- 
tions. One is, {(?q = d/8, g* = 0}, while the other set is 
given by, 



FIG. 1. The figure shows the flows in the coupling con- 
stants space for < p < d/2. The below one is for p > d/2. 



(10) 



which is stable in one cigendirection but unstable in 
the other eigendircction. 

The corresponding RG flow diagram is shown in Figure 
1. Therefore, for < p < d/2 the one-loop RG calcula- 
tion indicates that the system with uncorrelated disorder 
is unstable against long-range correlated disorder toward 
a new fixed point in the coupling constants space, for 
which there is a phase transition from delocalized to lo- 
calized states with increasing the disorder intensity. 

(ii) For p > d/2 there are two fixed points: the Gaus- 
sian fixed point which is stable on the go axis but unstable 
on the g p axis, and the 

non-trivial fixed point, {<?q = d/8, g* = 0}, which 
is unstable in all directions. The RG flow diagram for 
this case is shown in Figure 2. The implication is that, 
while the power- low correlated disorder is relevant, no 
new fixed point exists to one-loop order and, therefore, 
the long- wavelength behavior of the system is determined 
by the long-range component of the disorder. This means 
that for p > d/2 the waves are localized for any d. In 
addition, in both cases (i) and (ii) the system undergoes 
a disorder-induced transition when only the uncorrelated 
disorder is present. Let us mention that the above results 
are general so long as D p > (which is the only physically 
acceptable limit). For D p < the above phase space is 
valid for p > |(d + 1). 

To test these predictions, we have carried out numeri- 
cal simulations of the problem in both ID and 2D. Con- 
sider first the ID disordered systems. In this case, the 
waves are localized when the wave functions are of the 
form, ip(x) = f(x) cxp(— \x — ieo|/£)j where f(x) is a 
stochastic function which depends on the particular real- 
ization of the disordered chain, and £ is the localization 
length. Experimentally, the simplest ID model that ex- 
hibits wave localization is [18] a 15 m long steel wire with 
a 0.178 mm diameter, suspended vertically. The tension 
in the wire is maintained with a weight attached at its 



FIG. 2. The figure shows the flows in the coupling con- 
stants space for p > d/2. 




FIG. 3. The localization length £ for a disordered chain 
with 50,000 site and its dependence on the frequency ui 
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FIG. 4. Amplitude of the waves in 2D systems, for both 
randomly and uniformly-distributed and power-law correlated 
elastic constants. 



lower end. The function ij}{x, t) consists of transverse 
waves in the wire with an electromechanical actuator at 
one end of the wire. 

It was shown [18] that, even for very small deviations 
(less than 1%) from periodicity, the diagonal disorder 
(e.g., variations in the resonance frequencies of the oscil- 
lators) produces localization (which is in agreement with 
Furstenberg's theorem [19]), while variations (up to 13%) 
in the sizes of the masses (off-diagonal disorder) result in 
localization lengths that are much larger than the size 
of the system. This is in agreement with our theoretical 
prediction. 

To reproduce this result numerically and to calculate 
the localization length £, we used the transfer-matrix 
(TM) method [12]. Discrctizing Eq. (1) and writing 
down the result for site n of a linear chain yields, 

(u + \„)lp n + X n +l^n+2 - (A„ + l - A„)V>n+l = > (H) 

which can be rewritten in the recursive form 

(12) 

with 

— A ra + A n _2 A n _2 \ 

"x - ~y~ ) • (13) 

The localization length is then defined by, 

£(w) _1 = limTv^cx) AT - 1 1 -^jv / "00 1 > where N is the chain's 
length. For every realization of the disorder we computed 
ipN and, hence, £(u). We chose, ipo = V^i = l/V^j and 
averaged £ over a large ensemble of realizations for a fixed 
system size N and frequency lo. We then repeated this 
procedure for several values of N and lo. The extended 
states correspond to having, lirn/v^oo £/N = constant > 
1. As £ is also a function of lo, we chose, lo = 2ir\/Xo/N, 
which is the smallest mode of the system. Our RG anal- 
ysis indicates that this mode most likely passes through 
the ID disordered chain. 

The TM computations indicate that the coupling con- 
stant g(N) follows a finite-size scaling, g(N) = go + 1/N, 
where 170 is the coupling constant in the thermodynamic 
(N — > 00) limit; we find that, go ~ 0.117. In addition, 
when £ — > N, one can write, N tx (g — go)" 1 ' and obtain 
an estimate of the localization exponent v. Our analyti- 
cal results indicated [20] that in d dimensions, v = 1/d, 
which agrees with the TM calculations that yield v = 1. 
Figure 3 shows the results for the localization length £ as 
a function of lo. These results confirm the RG predictions 
for disordered linear chains. 

To further check the RG results, we also solved Eq. 
(1) in 2D using the finite-difference method with second- 
order discretization for both the space and time variables. 
Such approximations are acceptable as we work in the 
limit of low frequencies or long wavelengths. For short 



wavelengths we should use higher-order discretizations 
for the spatial variables [20]. A L x x L y grid was used 
with L x = 8000 and L y — 400. The parameter A(x), 
representing the local, gridblock-scale elastic constant, 
was distributed with a power-law correlation function, of 
the type considered above, with its spatial distribution 
generated using the midpoint displacement method [21]. 
We also simulated the case in which the local values of 
A(x) were uncorrelated and uniformly distributed, with 
the same variance as that of the power-law case. By in- 
serting a wave source on a line on one side of the grid, 
we simulated numerically 2D wave propagation through 
the system. Periodic boundary conditions were imposed 
in the lateral direction, which did not distort the nature 
of the wave propagation, as we used large system sizes. 
The decay in the amplitude of the wave is caused by scat- 
tering from heterogeneities of system generated by the 
distribution of the local elastic constants. The accuracy 
of the solution was checked by considering the stability 
criterion and the wavelength of the source [22] , and using 
higher-order finite-difference discretizations. To compute 
the amplitude decay in the medium, we collected the nu- 
merical results at 80 receivers (grid points), distributed 
evenly throughout the grid, along the direction of wave 
propagation. The results were averaged over 32 realiza- 
tions of the system. 

Figure 4 presents the decay in the wave amplitude 
through the uniformly random medium, and that of a 
medium with a nondecaying power-law correlation func- 
tion for the local elastic constants A(x), with p = 1.3 and 
1.8. The wave amplitudes for the correlated cases decline 
much faster than those in the uniformly random medium. 
In particular, for p = 1.3, which corresponds to negative 
correlations (that is, a large local elastic modulus is likely 
to be neighbor to a small one, and vice versa), the ampli- 
tude decreases rather sharply. These results confirm the 
RG predictions for 2D systems. 

In summary, we show that, depending on the nature 
of disorder, acoustic waves in strongly disordered media 
can be localized or delocalized in any dimensions. In 
particular, they can be extended in disordered ID sys- 
tems if the correlation function for the distribution of 
the local elastic constants is of nondecaying power-law 
type, and that the waves are localized in any dimension 
if the exponent p of the power-law correlation function 
is larger than d/2. These results, which contradict the 
generally-accepted view that off-diagonal disorder has a 
much weaker effect on localization than the diagonal dis- 
order, have important practical implications. For exam- 
ple, in order for seismic records to contain meaningful in- 
formation on the geology and content of a natural porous 
formation of linear size L, the localization length £ must 
be larger L. Otherwise, propagation and scattering of 
such waves can provide information on the formation only 
up to length scale £; one cannot obtain meaningful in- 
formation at larger length scales [23]. The localization 



length £ is, clearly, a function of the dimensionality of 
the system, the exponent p and amplitude D p , and other 
relevant physical parameters of the system. Its determi- 
nation remains a major task [20]. 
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